#ifndef XIBOV4_HPP
#define XIBOV4_HPP

#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <memory>
#include <omp.h>
#include <openblas/cblas.h>
#include <openblas/lapacke.h>

namespace XIBO
{
#ifndef XIBO_SCALAR
#define XIBO_SCALAR

#endif // XIBO_SCALAR

#ifndef XIBO_VECTOR
#define XIBO_VECTOR
    template <typename T>
    class Vector
    {
    private:
        int size_{};
        std::unique_ptr<T[]> data_{};

    public:
        Vector() = default;
        ~Vector()
        {
            if (data_ != nullptr)
                data_.reset();
        }
        explicit Vector(const int &n) : size_(n), data_(std::make_unique<T[]>(n)) {}
        Vector(const int &n, const T &a) : size_(n), data_(std::make_unique<T[]>(n))
        {
            std::fill_n(data_.get(), n, a);
        }
        Vector(const int &n, const T *a) : size_(n), data_(std::make_unique<T[]>(n))
        {
            std::copy_n(a, n, data_.get());
        }
        Vector(const Vector &rhs) : size_(rhs.size_), data_(std::make_unique<T[]>(rhs.size_))
        {
            std::copy_n(rhs.data_.get(), rhs.size_, data_.get());
        }
        Vector &operator=(const Vector &rhs)
        {
            if (this != &rhs)
            {
                Vector temp(rhs);
                std::swap(*this, temp);
            }
            return *this;
        }
        Vector(Vector &&rhs) noexcept
        {
            std::swap(*this, rhs);
        }
        Vector &operator=(Vector &&rhs) noexcept
        {
            if (this != &rhs)
            {
                std::swap(*this, rhs);
            }
            return *this;
        }
        T *begin()
        {
            return data_.get();
        }
        const T *begin() const
        {
            return data_.get();
        }
        int size() const
        {
            return size_;
        }
        T &operator()(const int &i)
        {
            if (i >= size_)
            {
                throw std::out_of_range("Index out of range");
            }
            return data_[i];
        }
        const T &operator()(const int &i) const
        {
            if (i >= size_)
            {
                throw std::out_of_range("Index out of range");
            }
            return data_[i];
        }
    };

#endif // XIBO_VECTOR

#ifndef XIBO_MATRIX
#define XIBO_MATRIX
    template <typename T>
    class Matrix
    {
    private:
        int row_{};
        int col_{};
        std::unique_ptr<T[]> data_{};

    public:
        Matrix() = default;
        ~Matrix()
        {
            if (data_ != nullptr)
            {
                data_.reset();
            }
        }
        explicit Matrix(const int &m, const int &n) : row_(m), col_(n), data_(std::make_unique<T[]>(m * n)) {}
        Matrix(const int &m, const int &n, const T &a) : row_(m), col_(n), data_(std::make_unique<T[]>(m * n))
        {
            std::fill_n(data_.get(), m * n, a);
        }
        Matrix(const int &m, const int &n, const T *a) : row_(m), col_(n), data_(std::make_unique<T[]>(m * n))
        {
            std::copy_n(a, m * n, data_.get());
        }
        Matrix(const Matrix &rhs) : row_(rhs.row_), col_(rhs.col_), data_(std::make_unique<T[]>(rhs.row_ * rhs.col_))
        {
            std::copy_n(rhs.data_.get(), rhs.size(), data_.get());
        }
        Matrix &operator=(const Matrix &rhs)
        {
            if (this != &rhs)
            {
                Matrix temp(rhs);
                std::swap(*this, temp);
            }
            return *this;
        }
        Matrix(Matrix &&rhs) noexcept
        {
            std::swap(*this, rhs);
        }
        Matrix &operator=(Matrix &&rhs) noexcept
        {
            if (this != &rhs)
            {
                std::swap(*this, rhs);
            }
            return *this;
        }
        T *begin()
        {
            return data_.get();
        }
        const T *begin() const
        {
            return data_.get();
        }
        int row() const
        {
            return row_;
        }
        int col() const
        {
            return col_;
        }
        int size() const
        {
            return row_ * col_;
        }
        T &operator()(const int &i, const int &j)
        {
            if (i >= row_ || j >= col_)
            {
                throw std::out_of_range("Index out of range");
            }
            return data_[i * col_ + j];
        }
        const T &operator()(const int &i, const int &j) const
        {
            if (i >= row_ || j >= col_)
            {
                throw std::out_of_range("Index out of range");
            }
            return data_[i * col_ + j];
        }
    };

#endif // XIBO_MATRIX

#ifndef XIBO_SOLVER
#define XIBO_SOLVER
    template <typename T>
    Vector<T> LapackeSolver(const Matrix<T> &A, const Vector<T> &b)
    {
        Matrix<T> A_(A);
        Vector<int> ipiv(A.row());
        Vector<T> x(b);
        LAPACKE_dgesv(LAPACK_ROW_MAJOR, A.row(), 1, A_.begin(), A.col(), ipiv.begin(), x.begin(), 1);

        return x;
    }

#endif // XIBO_SOLVER

#ifndef XIBO_FILEOPERATE
#define XIBO_FILEOPERATE

    template <typename T>
    void Print(const T &a)
    {
        std::cout << "Scalar: " << typeid(T).name() << std::endl
                  << a << std::endl;
    }

    template <typename T>
    void Print(const T *a, int size)
    {
        std::cout << "Array1D: " << typeid(T).name() << size << std::endl;
        std::cout << "[ ";
        for (int i = 0; i < size; ++i)
        {
            std::cout << *a++ << " ";
        }
        std::cout << "]" << std::endl;
    }

    template <typename T>
    void Print(const T *a, int row, int col)
    {
        std::cout << "Array2D: " << typeid(T).name() << row << "x" << col << std::endl;
        for (int i = 0; i < row; ++i)
        {
            std::cout << "[ ";
            for (int j = 0; j < col; ++j)
            {
                std::cout << *a++ << " ";
            }
            std::cout << "]" << std::endl;
        }
    }

    template <typename T>
    void Print(const Vector<T> &vec)
    {
        std::cout << "Vector: " << typeid(T).name() << vec.size() << std::endl;
        for (int i = 0; i < vec.size(); ++i)
        {
            std::cout << "[ " << vec(i) << " ]" << std::endl;
        }
    }

    template <typename T>
    void Print(const Matrix<T> &mat)
    {
        std::cout << "Matrix: " << typeid(T).name() << mat.row() << "x" << mat.col() << std::endl;
        for (int i = 0; i < mat.row(); ++i)
        {
            std::cout << "[ ";
            for (int j = 0; j < mat.col(); ++j)
            {
                std::cout << mat(i, j) << " ";
            }
            std::cout << "]" << std::endl;
        }
    }
#endif // XIBO_FILEOPERATE

#ifndef XIBO_NORM
#define XIBO_NORM
    template <typename T>
    double l2Norm(const Vector<T> &vec)
    {
        double result = 0;
        for (size_t i = 0; i < vec.size(); ++i)
        {
            result += std::pow(vec(i), 2);
        }
        return std::sqrt(result);
    }

    template <>
    double l2Norm<double>(const Vector<double> &vec)
    {
        return cblas_dnrm2(vec.size(), vec.begin(), 1);
    }

    template <>
    double l2Norm<float>(const Vector<float> &vec)
    {
        return cblas_snrm2(vec.size(), vec.begin(), 1);
    }
#endif // XIBO_NORM

}

#endif // XIBOV4_HPP